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Abstract 

A system of ODE's is used to attempt an approximation of the 
dynamics of two delayed coupled FitzHugh-Nagumo excitable units, 
described by delay-differential equations. It is shown that the codi- 
mension 2 generalized Hopf bifurcation acts as the organizing center 
for the dynamics of ODE's for small time-lags. Furthermore, this is 
used to explain important qualitative properties of the exact dynamics 
for small time-delays. 
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Introduction 



Dynamics of a pair of excitable systems with time-delayed coupling is 
quite different from the dynamics with the instantaneous coupling. Also, the 
behaviour of coupled excitable systems differs from that of coupled oscilla- 
tors with the same coupling. A prime example of the type II excitable be- 
haviour (see [Izhikevich, 2000]) has been the system introduced by FitzHugh 
[FitzHugh 1955] and Nagumo et all [Nagumo et all, 1962], as an approxima- 
tion of the Hodgine-Haxley model of the nerve cell membrane. One form of 
the FHN equations is (see [Murray, 1992]): 

i: = — + {a + l)x'^ — ax — y + I, 

y = hx-^y (1) 

As is well known the system (1) can operate in two regimes depending on 
the external current /, as an excitable system if / = and 

4- <(«-!)', (2) 
7 

with the stable stationary solution as the only attractor, or as a relaxation 
oscillator / = const ^ 0, with the stable limit cycle as the only attractor. 

Two delayed coupled FHN excitable systems with delayed coupling given 
by the following equations: 

Xi — —x\-\-{a-\-l)x\ — ax\ — yi-\-ciaj\'^{x2)-, 

yi = bxi - 7yi, 

X2 = —xl + {a + l)xl — ax2 — y2 + cta,n~^{x'[), 

y2 = bx2 - 7^2, (3) 

where x'^{t) — x{t — r), have been recently analyzed in [Buric&Todorovic, 
2003]. Motivation, background and the relevant hterature is discussed in 
detail in [Buric&Todorovic, 2003] and will not be repeated here. 

Besides different types of oscillations induced by coupling and changed by 
the time-delay, the system displays two different types of excitable behaviour. 
The first one is described by a single stationary solution at Eq = (0, 0, 0, 0) 
as the only attractor. The other type is described by two coexisting attrac- 
tors, the stable stationary solution Eq and a stable limit cycle correspond- 
ing to periodic exactly synchronous excitations of the two units. The later 
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regime occurs only in a specific, relatively small, domain in the parameters 
(c, r) plane. Buric&Todorovic obtained bifurcation curves in (c, r) plane by 
solving the characteristic equation of (3) for the non-hyperbolic roots. The 
bifurcation curves served as a guide for the analyzes of the system by ex- 
tensive numerical calculations, which resulted in the classification of possible 
excitable and oscillatory dynamics. 

In this letter we would like to report some analytic results aimed at bet- 
ter understanding of the types of bifurcations which are relevant for the two 
types of the excitable behaviour. To this end we shall analyze the bifurcations 
of the stationary solution that occur in a system related to the equations (3) 
in the following way. Both types of excitable behaviour happen for relatively 
small time-lags, and such a small time-lag induces the responsible bifurca- 
tions. Thus it might be justified to replace the time-delayed argument of the 
coupling function in (3) by the following approximation: 

f{xit-T))^f{x-Tx). (4) 

and approximate the delay-differential (DDE) by ordinary differential equa- 
tions. The approximate system is given by: 

Xi = —xl + (a + l)xl — axi — yi + 

+ ctan~^(a;2 — t{—xI + {a + l)xl — ax2 — I/2 + ctan~^(a;i))), 

yi = bxi - 71/1, 

X2 = —xl + {a + l)xl - ax2 - 1/2 + 

-|- ctan~^(a;i — t{—xI + {a + l)xl — axi — yi + ctan~^(a;2))), 

y2 = bx2 - 72/2, (5) 

Validity of the approximation (4) and (5) is not apriori justified even 
for small r. Such type of questions have been analyzed before in general 
and for specific examples (see [Meinardus & Nuenberg, 1985]). More recent 
example of such an analyzes, for a particular system, is given in the reference 
[Faro & Velasco, 1997], where the approximation has been investigated using 
a predator-prey equations with a time-delay by comparison of numerically 
obtained bifurcation curves. 

The main result of our analyzes is that the bifurcation which acts as 
the organizing center for the dynamics of the system (5) for small r is the 
codimension 2 generalized Hopf (Bautin) bifurcation. Furthermore, the bi- 
furcation occurs for quite a small time-lag where the bifurcation curves of the 
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exact and approximate system almost coincide and the dynamics is qualita- 
tively the same. This fact is used to explain the occurrence of the two types 
of excitable behaviour. 

The approximate and the exact systems will be abbreviated as X — 
T{X,X'^) and X — Tapp{X) where X e represents the collection of 
coordinates Xi,yi,X2,y2- 



Bifurcations of the stationary solution 



We shall restrict attention only on such values of the parameters a, b, 7, c 
that the system J^app for r = has only one stationary solution. This occurs 
if 

c < Ci = a + 6/7. (6) 

Futhermore, each unit is excitable when decoupled, i.e. the condition (2) is 
assumed satisfied. We fix the parameters a, b and 7 to some arbitrary such 
value, and consider the bifurcations of Eq that could occur as the parameters 
c and T are varied. The bifurcation set Beq is defined as the set of (c, r) 
values such that the stationary solution Eq is not hyperbolic. By an abuse of 
notation, we shall use the same symbol Beq for the part of Beq C R+ x R+ 
satisfying (6), i.e. 



Beo = {(c,t) e R+ X R+|ReA(c,T) = 0,c < ci} , 

where A is any root of the characteristic polynomial. 

The set Beq consists of two line segments in (c, r) plane. Namely: 



(7) 



B 



= {(c,r)|r = l/c,c< ci}|j|(c, 

= ^Eo;p[J^Eo;H- 

Indeed, the linear part of (5) 



T T 



c — a — ^ 
c{c — a) 



,c e (a,ci) 



(8) 
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where 



F — —a — c^T, D — c + car, E 



CT 



(9) 



(10) 
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implies the following characteristic equation: 

A(A) = Ai(A)A2(A) = 

where 



(11) 



Ai(A) = X'^ + {-f - F - D)X + b- F-f --fD -bE 
A2(A) = + (^- F- D)X + b- F-f + -fD + bE 



(12) 



with solutions 



Al,2 - 
^3,4 



-7 



+ F + D± ^(7 + F + L>)2 - 46(1 - £;)] /2, Ai(Ai,2) = 0, 



-7 + F - D ± a/(7 + F - - 46(1 + E) 



/2,A2(A3,4) = 0. (13) 



A nonhyperbolic root can be either equal to zero or pure imaginary. In the 
first case, 



Ai(0) = 44> b-F-f-gD-bE = 



(14) 



which defines the line segment Beq-p- The second factor of the characteristic 
polynomial has no zero roots for any positive c and r. In the second case 

Ai(w) = 0, i;>0 <^ -v^ + {^-F-D)iv + b-F-f--fD-bE^O 

^ = 6 - F7 - 7F) - 6F > and (7 - F - F))^; = 0. 

If c e (a, ci) then from the last condition we obtain the line segment Beq-h 

and V — ^''7(07 + 6 — crf)/{c — a). On the other hand, if c ^ (o, ci) there 
is no pure imaginary solution of Ai = 0. Furthermore A2 = has no pure 
imaginary solutions for any positive r and c. Thus, Beq is indeed given by 
(8). It is illustrated in figure la. The point where Beq intersects the c axis 
is denoted by cq. Thus: 

Co = a + 7. (15) 

The type of bifurcations occurring for the parameters (c, r) in Beo;p and 
Beq-h are described in the following two theorems. 

Theorem 1 The system J^app has a pitchfork bifurcation for any (c, r) e 

Beo;p- 
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Proof: The linear part A on Beq-^p has a simple zero eigenvalue. The 
type of bifurcation at r = 1/c is analyzed by reducing the system on the 
corresponding center manifold with the parameter e = t — 1/c. As we shall 
see, it is enough to consider the extended system, given by (5) and e = 
expended up to the third order in xi,X2,yi,y2, 

xi = -{a + c)xi - yi + {a + c)x2 + y2 + (a + '^){xj - xl) 

+ (c/3 — l){xl — X2) — axl — xly2 + cx\xi — c{a + l)exl + 

+ acex2 + cey2 — (?exi 

yi = hxi - 71/1 

X2 = {a + c)xi + yi - {a + c)x2 - y2 - {a + l){x\ - xi) 

— (c/3 — l)(a;'^ — X2) — ax\ — x\yi + cx\x2 — c{a + l)exl + 

+ acexi + ceyi — (?ex2 

y2 = bx2 - 71/2 

e = 0. (16) 

The center manifold with the parameter e of the system (16), in the new 
coordinates {x, y, z, t, e) related to the old ones by 

Xi — X — z — t, 

= bx/-f + y-bz/{-f + X3)-bt/{-f + X4), 

X2 = x + z + t, 

y2 = 6x/7 + y + 6^/(7 + A3) + 6V(7 + A4), 



A3,4 = (-2(a + c) -7±y'(2a + 2c-7)2-86)/2, (17) 

is: 

W^^iO) ^ {{x,y,z,t,e)\y = hi{x,e), z ^ h2{x,e),t ^ h3{x,e); 

hi{0, 0) = 0, ^^(O, 0) = 0, i = 1, 2, 3},(18) 

where 

hi(x,e) = --^(a + 6/7 - c)a;e + -^(a + 6/7 - c)a;^ + —(a + l)a;^e + ... , 
h2{x,e) = 0, 

hs{x,e) = 0. (19) 
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Restriction of (16) on the center manifold (18) is given by: 

X = F{x, e) = c(a + 6/7 — c)xe — {a + b/j — c)x'^ — c(a + l)x'^e + . . . , (20) 



and satisfies: 



dF{Q,Q) _^ d'F{Q,0) _^ 
de ' 

92^(0, 0) 



dedx 
d^F(0, 0) 



c{a + b/-f - c) 7^ 0, 
= -6(a + 6/7 - c) ^ 0. 



These are the sufficient and necessary conditions for the pitchfork bifurcation. 
The system (5), under the condition (6) is, in a neighborhood of {x, e) = (0, 0) 
locally topologically equivalent to x — ex — x^ (see [Arrowsmith, 1990], [ 
Kuznetsov, 1995]). Thus, if (6) is satisfied then for r ~ 1/c, and r < 1/c 
the stationary solution Eq is stable. For r > 1/c the stationary point Eq is 
unstable but there are two new stable stationary solutions. ^ 

Theorem 2 For the parameter values (c, r) e Beo;h the system J^^pp 
has either the supercritical Hopf or the subcritical Hopf or the generafized 
Hoph bifurcation. Furthermore, there are such values of a, b and 7 that the 
value cb for which the system has the generalized Hopf bifurcation satisfies 

Cb G (co,Ci). 

Proof: For the parameters in Beq,h the matrix A has a pair of purely 
imaginary eigenvalues Ai_2 = -^iv, v > Q and no other nonhyperbolic eigen- 
values. Furthermore, for (c, r) e Beo;h 

, rfReAi2, ld(-7 + F + L)), , , ,^ 

Thus, (c, r) G Beq;h corresponds to the Hopf bifurcation. The type of the 
Hopf bifurcation is determined by studying the normal form of the system on 
the two-dimensional invariant center manifold. To obtain the normal form 
we use the procedure introduced in [Coullet & Spiegel 1983], and applied 
by Kuznetsov [Kuznetsov, 1997] to obtain the relevant coefficients in normal 
forms of all codimension 1 and 2 bifurcations of stationary solutions of ODE. 

As we shall see, it is enough to start with the system J^app expanded up 
to the terms of the third order. 

1 1 

X = AX + -^^app,2{X^ X) -\- -J-app,'i{X^ X, X), 
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where 



/ (a + \)x\ — c{a + 1)txI \ 


(a + — c(a + 1)txI 




and 



( {c^t/S - l)xl + (cT - c/3 - caT)xl - CTx\yi + c^tx2x\ \ 



(?T I?) — + (cT — c/3 — carjxf — CTx\y\ + c^Ta;2a;i 





First introduce a complex eigenvector Q e of A, i.e. AQ — ivQ with 
components 



/ 1 \ 

(c — a)(l — it;/7) 



V(c-a)(l-iT;/7)/ 
and the corresponding eigenvector of A^: A^P — —ivP 



/ 



Av{c—a) 
v+i'y 
4v 

\ / 

normaUzed to < P,Q >= P^Q = 1. Vectors Q and Q (the complex conjugate 
of Q) form a basis in the center-subspace of A, so any vector R e E*^ can 
be written as i? = aQ + aQ where a —< P, R >E C^. The relation between 
the original system X = J-'app and the the complex normal form of the system 
on the center manifold X — H[a, a) of the following form: 



a 



iva + lia\a\^ + l2a\a\^ + 0{\a\^) 



(22) 



is contained in the corresponding homological equation: 

dH . dH ^ ^ 

-Q^OL + —a = Tapp\li\pL, OL)) 

Substituting the Taylor expansions of the transformation H and the system 
Tapp into the homological equation, and collecting the terms with the same 
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order one gets the coefficients Zi, Z2 ■ ■ ■ in tfie normal form. Tfie tliird order 
coefficient li is given by: 



h = ^Re < P, J'app,3(Q,<5,<5) +^app,2((5, (2iv/4 - ^) ^:^app,2(Q,(3)) 



(c + 3)- 



7 



+ c - a - 7 + 



27'(a + l) 



wfiere 



c — a [c — a){b + a-y — 07) 

7c^ + E(a, b, 7)c^ + C(a, fe, 7)0 + /^(a, 6, 7) 
2{a — c){b + aj — cj) ' 



B{a,b,-f) = -fc-3a7-272 

C(a,6,7) = 2a6 + 3a^7 + 267 - 37^ + 807^ 

D(a, 6, 7) = —a^b — c?^ + 867 — 067 — 27^ — 07^ — 3a'^7'^. 



(28) 



(24) 



If tlie tliird order coefficient /i 7^ the system is locally smoothly orbitally 
equivalent to the system 

a — iva + Ziq;|q;|^. 

Since everywhere on Beq,h we have < , the values of the parameters 
(c, r) G Beo,h such that h < 0, imply the supercitical Hopf bifurcation, 
if Zi > the bifurcation is subcritical, and if Zi = the bifurcation is of 
the generalized Hopf type. The denominator of li is always negative in the 
considered interval c G (a, Ci = a + b/'~f). The numerator of li always has 
at least one real zero cb- Thus, there is at least one value of c = such 
that the bifurcation is of the generalized Hopf type. The following choice 
of a, b and 7: a = 0.25, b — ^ — 0.02 is an example of the values such that 
Cb = 0.289024 G (co,ci) = (0.27,1,27). In this case all three alternatives 
occur as c is varied in the interval (co, Ci). 

In the case cb G (cq, Ci) all three types of Hopf bifurcation occur. In 
general, the line segment of codimension 1 subcritical Hopf bifurcations, joins 
with the line segment of codimension 1 supercritical Hopf bifurcations at 
some point (c^. tb) of the codimension 2 generalized Hopf bifurcation. It is 
important to point out that in all numerical test that we have performed for 
the values of a, b and 7 such that the isolated unit shows excitable behaviour. 
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we have observed precisely the situation where cb £ (cq, Ci) and actually is 
quite close to Cq. 

As is well known from the theory of the codimension 2 generalized Hopf 
bifurcation, besides the two lines of Hopf bifurcations of the stationary solu- 
tion, there is also a line of fold limit cycle bifurcations emanating from the 
point {cb,tb)- For the parameter values in between the lines of subcritical 
Hopf and fold limit cycle bifurcations the system Tapp has a stable station- 
ary solution surrounded by an unstable limit cycle, which is surrounded by 
a stable limit cycle, all three lying in the manifold xi = X2,yi = 2/2- The 
generalized Hopf bifurcation at {cb,tb) is illustrated in figure lb. 

Finally, let us remark that the theorems 1 and 2 are probably correct if 
the tan~^ coupling function is replaced by any function of the sigmoid form, 
although we do not gave the proof in the general case because the algebra 
gets rather tedious. The only conditions that should be required are that: 
/(O) = 0, /' (0) > 0, r (0) = and (0) ^ 0. 

Approximate vs exact system 

The bifurcation set in (c, r) plane of the exact system under the same 
conditions (2) and (6) on the parameters a, b, 7 and c, was obtained before in 
[Buric & Todorovic, 2003]. Using the bifurcation set and the numerical test, it 
was conjectured that there is a domain in (c, r) that corresponds to the death 
of oscillations due to time-delay. On the bases of the numerical evidence it 
was also conjectured that the bifurcation mechanism beyond the oscillator 
death is more complicated than commonly found in oscillators coupled by 
diffusion with delay, i.e. the inverse supercritical Hopf bifurcation (see [Reddy 
et. all., 1999]). 

We would like to use the results about the bifurcations of the approximate 
system J^^pp in order to discuss the time-delay death of oscillations that 
have been induced by coupling of excitable systems. To this end, we first 
compare the bifurcation sets in (c, r) plane of the exact end the approximate 
system. The two sets arc illustrated in figure 2. The curves denoted Ti and 
T2 correspond to the first and second factor of the characteristic equation of 
the exact system, and there dReXi^2/dT < on ti and dReX3^4^/dT < on 
T2, where Ai^2,3,4 are the solutions of the exact characteristic equation. 

Few observations are in order. Firstly, the approximate systems badly 
fails to describe the bifurcations of the exact system for most values of the 
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time-lag. Qualitative agreement between the dynamics of the two systems 
is obtained only for very small values of the time-lag. The entire family of 
bifurcations due to nonhyperbolicity of the second factor of the characteris- 
tic equation of the exact system is missing. The smallest time-lag for these 
bifurcations to occur, is to large to be capture by the approximation, for all 
values of c < Ci. Furthermore, the approximate system has a line of pitch- 
fork bifurcations, destabilizing Eq and introducing two new stable stationary 
solutions for large r and any c. There is no analogous situation in the exact 
system. 

On the other hand, there is a small but important domain of (c, r) where 
the bifurcation curves of the exact system are well approximated by the ap- 
proximate one. This is the domain precisely around the generalized Hopf 
bifurcation, and is indicated in the figure 2. Thus, we can use theorem 2, 
about the bifurcations of the approximate system, to explain the dynamics 
of the exact system near the bifurcation line ri(c). This also supports the 
conjecture that the mechanism involved in the death of oscillations due to 
time-delay in the exact system involves the line of subcritical Hopf and the 
line of fold limit cycle bifurcations organized by the generalized Hopf bifur- 
cation. Consider the system for c > Cb and zero or small time-lag, when it 
consists of two exactly synchronized oscillators (compare figures 3 and 4). 
The only attractor is the limit cycle in Xi — X2, yi — y2 plane. Increasing 
the time-lag leads to the subcritical Hopf bifurcation for (c, r) e Beq,h when 
the stationary solution becomes stable and an unstable limit cycle is created 
in the same plane as the stable limit cycle. The system is bi-stable with the 
stable stationary state and periodic excitations described by the stable limit 
cycle. The unstable hmit cycle acts as a threshold. Further increase of the 
time lag leads to the disappearance of the two limit cycles in the fold limit 
cycle bifurcation. This corresponds to the death of oscillatory excitations. 
Upon further increase of r the approximate system hits the line of pitch- 
fork bifurcations, the stationary point becomes unstable and there are two 
new stable stationary solutions. Such dynamics does not occur in the exact 
system for any value of the time-lag. The sequence of different attractors 
obtained for fixed c and successively larger r for the approximate system, 
illustrated in figure 3, corresponds qualitatively to the sequence in the exact 
system illustrated in figure 4. Large time-lags, introduce qualitatively differ- 
ent dynamics (figures 3d and 4d). The qualitative correspondence between 
the exact system and the approximation is lost. 
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Summary and conclusions 

We have performed an analyzes of bifurcations of the stationary solution 
of a model of two coupled FitzHugh-Nagumo excitable systems. The model 
ODE's originate as a small time-lag approximation of the DDE's which ex- 
plicitly include the time-delay in the transmission of excitations between the 
units. 

The main results are given in the two theorems 1 and 2. The second 
theorem identifies the degenerate Hopf bifurcation as the main organizing 
center that is enough to explain all qualitative features of the dynamics for 
small time-lags. There are three possible types of dynamics in this domain. 
The system could be excitable, with the stable stationary solution as the 
only attractor, or oscillatory, when the limit cycle is the only attractor, or, 
finally, the system could be bi-stable. In the last case there are the stable 
stationary solution and the stable hmit cycle. These are the three types of 
the dynamics that have been observed also in the exact system of DDE for 
small time-lags. The first theorem determines the boundary in the (c, r) 
plain beyond which the dynamics of the approximate system is qualitatively 
different from anything that occurs in the exact system for the considered 
parameter values. 

Our analyzes is carried out using an explicit coupling function. However, 
the results would probably be the same for any coupling of the same form with 
the function f\x) satisfying /(O) = 0, /'(O) > 0, /"(O) = and /"'(O) ^ 0. 
On the other hand, different type of coupling, for example of the diffusive 
form, would result in different bifurcations and dynamics. 

There is yet another set of codimension 2 bifurcations that occur in the 
exact system for larger time-delays. They happen at the intersection of 
Hopf bifurcation curves, and are quite important for the properties of the 
oscillatory dynamics at larger time-lags. This Hopf-Hopf bifurcations are 
not captured by the finite dimensional approximation by ODE's. In order to 
analyze them an infinite dimensional generalization of the method applied in 
theorem 2 (see [Faria & Magelhes, 1995], [Schayer & Campbell, 1998]) could 
be applied. 

Acknowledgements This work is supported by the Serbian Ministry of 
Science contract No. 101443. 
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FIGURE CAPTIONS 

Figure 1: (a) Bifurcation set of the approximate system, the encircled 
area is enlarged in (b); (b) Dynamics near the codimension 2 generalized 
Hopf bifurcation. 

Figure 2: Bifurcation sets of the exact (thick lines) and the approximate 
(thin lines) systems. Ti and T2 are Hopf bifurcation curves of the exact 
system, and TH^app and Tp^app ^-re Hopf and pitchfork bifurcation lines of the 
approximate system. 

Figure 3: Projections on {xi, X2) of typical orbits approaching the possible 
attractors of the approximate system: a) One stable limit cycle (symmetric), 
b) stable stationary solution and the stable limit cycle (symmetric), c)one 
stable stationary solution, d) two stable stationary solutions. 

Figure 4: Projections on {xi,X2) of typical orbits approaching the possible 
attractors of the exact system, a) One stable limit cycle (symmetric), b) 
stable stationary solution and the stable limit cycle (symmetric), c)one stable 
stationary solution, d) one stable (asymmetric) limit cycle. 
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